home *** CD-ROM | disk | FTP | other *** search
/ Aminet 1 (Walnut Creek) / Aminet - June 1993 [Walnut Creek].iso / usenet / sources / volume2 / aplictns / matlab / src.3 < prev    next >
Internet Message Format  |  1988-11-02  |  49KB
  1. Path: xanth!nic.MR.NET!hal!cwjcc!mailrus!ulowell!page
  2. From: page@swan.ulowell.edu (Bob Page)
  3. Newsgroups: comp.sources.amiga
  4. Subject: v02i043:  matlab - matrix laboratory, Part03/11
  5. Message-ID: <10018@swan.ulowell.edu>
  6. Date: 2 Nov 88 21:41:59 GMT
  7. Organization: University of Lowell, Computer Science Dept.
  8. Lines: 1220
  9. Approved: page@swan.ulowell.edu
  10.  
  11. Submitted-by: strovink%galaxy-43@afit-ab.arpa (Mark A. Strovink)
  12. Posting-number: Volume 2, Issue 43
  13. Archive-name: applications/matlab/src.3
  14.  
  15. #    This is a shell archive.
  16. #    Remove everything above and including the cut line.
  17. #    Then run the rest of the file through sh.
  18. #----cut here-----cut here-----cut here-----cut here----#
  19. #!/bin/sh
  20. # shar:    Shell Archiver
  21. #    Run the following text with /bin/sh to create:
  22. #    src-3
  23. # This archive created: Wed Nov  2 16:20:21 1988
  24. cat << \SHAR_EOF > src-3
  25.    44 MRAT = IDINT(STKR(L))               
  26.       LRAT = IDINT(STKR(L-1))             
  27.       TOP = TOP - 1    
  28.       MSTK(TOP) = 0    
  29.       GO TO 99         
  30. C                      
  31. C     CHAR             
  32.    50 K = IABS(IDINT(STKR(L)))            
  33.       IF (M*N.NE.1 .OR. K.GE.ALFL) CALL ERROR(36)            
  34.       IF (ERR .GT. 0) RETURN              
  35.       CH = ALFA(K+1)   
  36.       IF (STKR(L) .LT. 0.0D0) CH = ALFB(K+1)                 
  37.       WRITE(WTE,51) CH                    
  38.    51 FORMAT(1X,'REPLACE CHARACTER ',A1)  
  39.       READ(RTE,52) CH                     
  40.    52 FORMAT(A1)       
  41.       IF (STKR(L) .GE. 0.0D0) ALFA(K+1) = CH                 
  42.       IF (STKR(L) .LT. 0.0D0) ALFB(K+1) = CH                 
  43.       MSTK(TOP) = 0    
  44.       GO TO 99         
  45. C                      
  46. C     DISP             
  47.    60 WRITE(WTE,61)    
  48.       IF (WIO .NE. 0) WRITE(WIO,61)       
  49.    61 FORMAT(1X,80A1)                     
  50.       IF (RHS .EQ. 2) GO TO 65            
  51.       MN = M*N         
  52.       TEXT = .TRUE.    
  53.       DO 62 I = 1, MN                     
  54.         LS = L+I-1     
  55.         CH = IDINT(STKR(LS))              
  56.         TEXT = TEXT .AND. (CH.GE.0) .AND. (CH.LT.ALFL)       
  57.         TEXT = TEXT .AND. (DFLOAT(CH).EQ.STKR(LS))           
  58.    62 CONTINUE         
  59.       DO 64 I = 1, M   
  60.       DO 63 J = 1, N   
  61.         LS = L+I-1+(J-1)*M                
  62.         IF (STKR(LS) .EQ. 0.0D0) CH = BLANK                  
  63.         IF (STKR(LS) .GT. 0.0D0) CH = PLUS                   
  64.         IF (STKR(LS) .LT. 0.0D0) CH = MINUS                  
  65.         IF (TEXT) CH = IDINT(STKR(LS))    
  66.         BUF(J) = ALFA(CH+1)               
  67.    63 CONTINUE         
  68.       WRITE(WTE,61) (BUF(J),J=1,N)        
  69.       IF (WIO .NE. 0) WRITE(WIO,61) (BUF(J),J=1,N)           
  70.    64 CONTINUE         
  71.       MSTK(TOP) = 0    
  72.       GO TO 99         
  73. C                      
  74. C     BASE             
  75.    65 IF (RHS .NE. 2) CALL ERROR(39)      
  76.       IF (STKR(L) .LE. 1.0D0) CALL ERROR(36)                 
  77.       IF (ERR .GT. 0) RETURN              
  78.       B = STKR(L)      
  79.       L2 = L           
  80.       TOP = TOP-1      
  81.       RHS = 1          
  82.       L = LSTK(TOP)    
  83.       M = MSTK(TOP)*NSTK(TOP)             
  84.       EPS = STKR(VSIZE-4)                 
  85.       DO 66 I = 1, M   
  86.          LS = L2+(I-1)*N                  
  87.          LL = L+I-1    
  88.          CALL BASE(STKR(LL),B,EPS,STKR(LS),N)                
  89.    66 CONTINUE         
  90.       CALL RSET(M*N,0.0D0,STKI(L2),1)     
  91.       CALL WCOPY(M*N,STKR(L2),STKI(L2),1,STKR(L),STKI(L),1)  
  92.       MSTK(TOP) = N    
  93.       NSTK(TOP) = M    
  94.       CALL STACK1(QUOTE)                  
  95.       IF (FIN .EQ. 6) GO TO 60            
  96.       GO TO 99         
  97. C                      
  98. C     LINES            
  99.    70 LCT(2) = IDINT(STKR(L))             
  100.       MSTK(TOP) = 0    
  101.       GO TO 99         
  102. C                      
  103. C     PLOT             
  104.    80 IF (RHS .GE. 2) GO TO 82            
  105.       N = M*N          
  106.       DO 81 I = 1, N   
  107.          LL = L+I-1    
  108.          STKI(LL) = DFLOAT(I)             
  109.    81 CONTINUE         
  110.       CALL PLOT(WTE,STKI(L),STKR(L),N,T,0,BUF)               
  111.       IF (WIO .NE. 0) CALL PLOT(WIO,STKI(L),STKR(L),N,T,0,BUF)                  
  112.       MSTK(TOP) = 0    
  113.       GO TO 99         
  114.    82 IF (RHS .EQ. 2) K = 0               
  115.       IF (RHS .EQ. 3) K = M*N             
  116.       IF (RHS .GT. 3) K = RHS - 2         
  117.       TOP = TOP - (RHS - 1)               
  118.       N = MSTK(TOP)*NSTK(TOP)             
  119.       IF (MSTK(TOP+1)*NSTK(TOP+1) .NE. N) CALL ERROR(5)      
  120.       IF (ERR .GT. 0) RETURN              
  121.       LX = LSTK(TOP)   
  122.       LY = LSTK(TOP+1)                    
  123.       IF (RHS .GT. 3) L = LSTK(TOP+2)     
  124.       CALL PLOT(WTE,STKR(LX),STKR(LY),N,STKR(L),K,BUF)       
  125.       IF (WIO .NE. 0) CALL PLOT(WIO,STKR(LX),STKR(LY),N,STKR(L),K,BUF)          
  126.       MSTK(TOP) = 0    
  127.       GO TO 99         
  128. C                      
  129. C     DEBUG            
  130.    95 DDT = IDINT(STKR(L))                
  131.       WRITE(WTE,96) DDT                   
  132.    96 FORMAT(1X,'DEBUG ',I4)              
  133.       MSTK(TOP) = 0    
  134.       GO TO 99         
  135. C                      
  136.    99 RETURN           
  137.       END
  138.               
  139.       SUBROUTINE MATFN6                   
  140. C                      
  141. C     EVALUATE UTILITY FUNCTIONS          
  142. C                      
  143.       DOUBLE PRECISION STKR(5005),STKI(5005)                 
  144.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP        
  145.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),HIO,RIO,WIO,RTE,WTE,FE           
  146.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)            
  147.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP          
  148.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,HIO,RIO,WIO,RTE,WTE,FE                  
  149.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN  
  150.       INTEGER SEMI,ID(4),UNIFOR(4),NORMAL(4),SEED(4)         
  151.       DOUBLE PRECISION EPS0,EPS,S,SR,SI,T                    
  152.       DOUBLE PRECISION FLOP,URAND         
  153.       LOGICAL EQID     
  154.       DATA SEMI/39/    
  155.       DATA UNIFOR/30,23,18,15/,NORMAL/23,24,27,22/,SEED/28,14,14,13/            
  156. C                      
  157.       IF (DDT .EQ. 1) WRITE(WTE,100) FIN  
  158.   100 FORMAT(1X,'MATFN6',I4)              
  159. C     FUNCTIONS/FIN    
  160. C     MAGI DIAG SUM  PROD USER EYE  RAND ONES CHOP SIZE KRON  TRIL TRIU         
  161. C       1    2    3    4    5    6    7    8    9   10  11-13  14   15          
  162.       L = LSTK(TOP)    
  163.       M = MSTK(TOP)    
  164.       N = NSTK(TOP)    
  165.       GO TO (75,80,65,67,70,90,90,90,60,77,50,50,50,80,80),FIN                  
  166. C                      
  167. C     KRONECKER PRODUCT                   
  168.    50 IF (RHS .NE. 2) CALL ERROR(39)      
  169.       IF (ERR .GT. 0) RETURN              
  170.       TOP = TOP - 1    
  171.       L = LSTK(TOP)    
  172.       MA = MSTK(TOP)   
  173.       NA = NSTK(TOP)   
  174.       LA = L + MAX0(M*N*MA*NA,M*N+MA*NA)  
  175.       LB = LA + MA*NA                     
  176.       ERR = LB + M*N - LSTK(BOT)          
  177.       IF (ERR .GT. 0) CALL ERROR(17)      
  178.       IF (ERR .GT. 0) RETURN              
  179. C     MOVE A AND B ABOVE RESULT           
  180.       CALL WCOPY(MA*NA+M*N,STKR(L),STKI(L),1,STKR(LA),STKI(LA),1)               
  181.       DO 54 JA = 1, NA                    
  182.         DO 53 J = 1, N                    
  183.           LJ = LB + (J-1)*M               
  184.           DO 52 IA = 1, MA                
  185. C           GET J-TH COLUMN OF B          
  186.             CALL WCOPY(M,STKR(LJ),STKI(LJ),1,STKR(L),STKI(L),1)                 
  187. C           ADDRESS OF A(IA,JA)           
  188.             LS = LA + IA-1 + (JA-1)*MA    
  189.             DO 51 I = 1, M                
  190. C             A(IA,JA) OP B(I,J)          
  191.               IF (FIN .EQ. 11) CALL WMUL(STKR(LS),STKI(LS),  
  192.      $           STKR(L),STKI(L),STKR(L),STKI(L))            
  193.               IF (FIN .EQ. 12) CALL WDIV(STKR(LS),STKI(LS),  
  194.      $           STKR(L),STKI(L),STKR(L),STKI(L))            
  195.               IF (FIN .EQ. 13) CALL WDIV(STKR(L),STKI(L),    
  196.      $           STKR(LS),STKI(LS),STKR(L),STKI(L))          
  197.               IF (ERR .GT. 0) RETURN      
  198.               L = L + 1                   
  199.    51       CONTINUE   
  200.    52     CONTINUE     
  201.    53   CONTINUE       
  202.    54 CONTINUE         
  203.       MSTK(TOP) = M*MA                    
  204.       NSTK(TOP) = N*NA                    
  205.       GO TO 99         
  206. C                      
  207. C     CHOP             
  208.    60 EPS0 = 1.0D0     
  209.    61 EPS0 = EPS0/2.0D0                   
  210.       T = FLOP(1.0D0 + EPS0)              
  211.       IF (T .GT. 1.0D0) GO TO 61          
  212.       EPS0 = 2.0D0*EPS0                   
  213.       FLP(2) = IDINT(STKR(L))             
  214.       IF (SYM .NE. SEMI) WRITE(WTE,62) FLP(2)                
  215.    62 FORMAT(/1X,'CHOP ',I2,' PLACES.')   
  216.       EPS = 1.0D0      
  217.    63 EPS = EPS/2.0D0                     
  218.       T = FLOP(1.0D0 + EPS)               
  219.       IF (T .GT. 1.0D0) GO TO 63          
  220.       EPS = 2.0D0*EPS                     
  221.       T = STKR(VSIZE-4)                   
  222.       IF (T.LT.EPS .OR. T.EQ.EPS0) STKR(VSIZE-4) = EPS       
  223.       MSTK(TOP) = 0    
  224.       GO TO 99         
  225. C                      
  226. C     SUM              
  227.    65 SR = 0.0D0       
  228.       SI = 0.0D0       
  229.       MN = M*N         
  230.       DO 66 I = 1, MN                     
  231.          LS = L+I-1    
  232.          SR = FLOP(SR+STKR(LS))           
  233.          SI = FLOP(SI+STKI(LS))           
  234.    66 CONTINUE         
  235.       GO TO 69         
  236. C                      
  237. C     PROD             
  238.    67 SR = 1.0D0       
  239.       SI = 0.0D0       
  240.       MN = M*N         
  241.       DO 68 I = 1, MN                     
  242.          LS = L+I-1    
  243.          CALL WMUL(STKR(LS),STKI(LS),SR,SI,SR,SI)            
  244.    68 CONTINUE         
  245.    69 STKR(L) = SR     
  246.       STKI(L) = SI     
  247.       MSTK(TOP) = 1    
  248.       NSTK(TOP) = 1    
  249.       GO TO 99         
  250. C                      
  251. C     USER             
  252.    70 S = 0.0D0        
  253.       T = 0.0D0        
  254.       IF (RHS .LT. 2) GO TO 72            
  255.       IF (RHS .LT. 3) GO TO 71            
  256.       T = STKR(L)      
  257.       TOP = TOP-1      
  258.       L = LSTK(TOP)    
  259.       M = MSTK(TOP)    
  260.       N = NSTK(TOP)    
  261.    71 S = STKR(L)      
  262.       TOP = TOP-1      
  263.       L = LSTK(TOP)    
  264.       M = MSTK(TOP)    
  265.       N = NSTK(TOP)    
  266.    72 CALL USER(STKR(L),M,N,S,T)          
  267.       CALL RSET(M*N,0.0D0,STKI(L),1)      
  268.       MSTK(TOP) = M    
  269.       NSTK(TOP) = N    
  270.       GO TO 99         
  271. C                      
  272. C     MAGIC            
  273.    75 N = MAX0(IDINT(STKR(L)),0)          
  274.       IF (N .EQ. 2) N = 0                 
  275.       IF (N .GT. 0) CALL MAGIC(STKR(L),N,N)                  
  276.       CALL RSET(N*N,0.0D0,STKI(L),1)      
  277.       MSTK(TOP) = N    
  278.       NSTK(TOP) = N    
  279.       GO TO 99         
  280. C                      
  281. C     SIZE             
  282.    77 STKR(L) = M      
  283.       STKR(L+1) = N    
  284.       STKI(L) = 0.0D0                     
  285.       STKI(L+1) = 0.0D0                   
  286.       MSTK(TOP) = 1    
  287.       NSTK(TOP) = 2    
  288.       IF (LHS .EQ. 1) GO TO 99            
  289.       NSTK(TOP) = 1    
  290.       TOP = TOP + 1    
  291.       LSTK(TOP) = L+1                     
  292.       MSTK(TOP) = 1    
  293.       NSTK(TOP) = 1    
  294.       GO TO 99         
  295. C                      
  296. C     DIAG, TRIU, TRIL                    
  297.    80 K = 0            
  298.       IF (RHS .NE. 2) GO TO 81            
  299.          K = IDINT(STKR(L))               
  300.          TOP = TOP-1   
  301.          L = LSTK(TOP)                    
  302.          M = MSTK(TOP)                    
  303.          N = NSTK(TOP)                    
  304.    81 IF (FIN .GE. 14) GO TO 85           
  305.       IF (M .EQ. 1 .OR. N .EQ. 1) GO TO 83                   
  306.       IF (K.GE.0) MN=MIN0(M,N-K)          
  307.       IF (K.LT.0) MN=MIN0(M+K,N)          
  308.       MSTK(TOP) = MAX0(MN,0)              
  309.       NSTK(TOP) = 1    
  310.       IF (MN .LE. 0) GO TO 99             
  311.       DO 82 I = 1, MN                     
  312.          IF (K.GE.0) LS = L+(I-1)+(I+K-1)*M                  
  313.          IF (K.LT.0) LS = L+(I-K-1)+(I-1)*M                  
  314.          LL = L+I-1    
  315.          STKR(LL) = STKR(LS)              
  316.          STKI(LL) = STKI(LS)              
  317.    82 CONTINUE         
  318.       GO TO 99         
  319.    83 N = MAX0(M,N)+IABS(K)               
  320.       ERR = L+N*N - LSTK(BOT)             
  321.       IF (ERR .GT. 0) CALL ERROR(17)      
  322.       IF (ERR .GT. 0) RETURN              
  323.       MSTK(TOP) = N    
  324.       NSTK(TOP) = N    
  325.       DO 84 JB = 1, N                     
  326.       DO 84 IB = 1, N                     
  327.          J = N+1-JB    
  328.          I = N+1-IB    
  329.          SR = 0.0D0    
  330.          SI = 0.0D0    
  331.          IF (K.GE.0) LS = L+I-1           
  332.          IF (K.LT.0) LS = L+J-1           
  333.          LL = L+I-1+(J-1)*N               
  334.          IF (J-I .EQ. K) SR = STKR(LS)    
  335.          IF (J-I .EQ. K) SI = STKI(LS)    
  336.          STKR(LL) = SR                    
  337.          STKI(LL) = SI                    
  338.    84 CONTINUE         
  339.       GO TO 99         
  340. C                      
  341. C     TRIL, TRIU       
  342.    85 DO 87 J = 1, N   
  343.          LD = L + J - K - 1 + (J-1)*M     
  344.          IF (FIN .EQ. 14) LL = J - K - 1  
  345.          IF (FIN .EQ. 14) LS = LD - LL    
  346.          IF (FIN .EQ. 15) LL = M - J + K  
  347.          IF (FIN .EQ. 15) LS = LD + 1     
  348.          IF (LL .GT. 0) CALL WSET(LL,0.0D0,0.0D0,STKR(LS),STKI(LS),1)           
  349.    87 CONTINUE         
  350.       GO TO 99         
  351. C                      
  352. C     EYE, RAND, ONES                     
  353.    90 IF (M.GT.1 .OR. RHS.EQ.0) GO TO 94  
  354.       IF (RHS .NE. 2) GO TO 91            
  355.         NN = IDINT(STKR(L))               
  356.         TOP = TOP-1    
  357.         L = LSTK(TOP)                     
  358.         N = NSTK(TOP)                     
  359.    91 IF (FIN.NE.7 .OR. N.LT.4) GO TO 93  
  360.       DO 92 I = 1, 4   
  361.         LS = L+I-1     
  362.         ID(I) = IDINT(STKR(LS))           
  363.    92 CONTINUE         
  364.       IF (EQID(ID,UNIFOR).OR.EQID(ID,NORMAL)) GO TO 97       
  365.       IF (EQID(ID,SEED)) GO TO 98         
  366.    93 IF (N .GT. 1) GO TO 94              
  367.       M = MAX0(IDINT(STKR(L)),0)          
  368.       IF (RHS .EQ. 2) N = MAX0(NN,0)      
  369.       IF (RHS .NE. 2) N = M               
  370.       ERR = L+M*N - LSTK(BOT)             
  371.       IF (ERR .GT. 0) CALL ERROR(17)      
  372.       IF (ERR .GT. 0) RETURN              
  373.       MSTK(TOP) = M    
  374.       NSTK(TOP) = N    
  375.       IF (M*N .EQ. 0) GO TO 99            
  376.    94 DO 96 J = 1, N   
  377.       DO 96 I = 1, M   
  378.         LL = L+I-1+(J-1)*M                
  379.         STKR(LL) = 0.0D0                  
  380.         STKI(LL) = 0.0D0                  
  381.         IF (I.EQ.J .OR. FIN.EQ.8) STKR(LL) = 1.0D0           
  382.         IF (FIN.EQ.7 .AND. RAN(2).EQ.0) STKR(LL) = FLOP(URAND(RAN(1)))          
  383.         IF (FIN.NE.7 .OR. RAN(2).EQ.0) GO TO 96              
  384.    95      SR = 2.0D0*URAND(RAN(1))-1.0D0                    
  385.            SI = 2.0D0*URAND(RAN(1))-1.0D0                    
  386.            T = SR*SR + SI*SI              
  387.            IF (T .GT. 1.0D0) GO TO 95     
  388.         STKR(LL) = FLOP(SR*DSQRT(-2.0D0*DLOG(T)/T))          
  389.    96 CONTINUE         
  390.       GO TO 99         
  391. C                      
  392. C     SWITCH UNIFORM AND NORMAL           
  393.    97 RAN(2) = ID(1) - UNIFOR(1)          
  394.       MSTK(TOP) = 0    
  395.       GO TO 99         
  396. C                      
  397. C     SEED             
  398.    98 IF (RHS .EQ. 2) RAN(1) = NN         
  399.       STKR(L) = RAN(1)                    
  400.       MSTK(TOP) = 1    
  401.       IF (RHS .EQ. 2) MSTK(TOP) = 0       
  402.       NSTK(TOP) = 1    
  403.       GO TO 99         
  404. C                      
  405.    99 RETURN           
  406.       END
  407.       SUBROUTINE MATLAB(INIT)             
  408. C     INIT = 0 FOR ORDINARY FIRST ENTRY   
  409. C          = POSITIVE FOR SUBSEQUENT ENTRIES                 
  410. C          = NEGATIVE FOR SILENT INITIALIZATION (SEE MATZ)   
  411. C                      
  412.       DOUBLE PRECISION STKR(5005),STKI(5005)                 
  413.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP        
  414.       INTEGER ALFA(52),ALFB(52),ALFL,CASE                    
  415.       INTEGER IDS(4,32),PSTK(32),RSTK(32),PSIZE,PT,PTZ       
  416.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),HIO,RIO,WIO,RTE,WTE,FE           
  417.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)            
  418.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP          
  419.       COMMON /ALFS/ ALFA,ALFB,ALFL,CASE   
  420.       COMMON /RECU/ IDS,PSTK,RSTK,PSIZE,PT,PTZ               
  421.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,HIO,RIO,WIO,RTE,WTE,FE                  
  422.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN  
  423. C                      
  424.       DOUBLE PRECISION S,T                
  425.       INTEGER EPS(4),FLOPS(4),EYE(4),RAND(4)                 
  426. C                      
  427. C     CHARACTER SET    
  428. C            0       10       20       30       40       50  
  429. C                      
  430. C     0      0        A        K        U   COLON  :  LESS   <                  
  431. C     1      1        B        L        V   PLUS   +  GREAT  >                  
  432. C     2      2        C        M        W   MINUS  -         
  433. C     3      3        D        N        X   STAR   *         
  434. C     4      4        E        O        Y   SLASH  /         
  435. C     5      5        F        P        Z   BSLASH \         
  436. C     6      6        G        Q  BLANK     EQUAL  =         
  437. C     7      7        H        R  LPAREN (  DOT    .         
  438. C     8      8        I        S  RPAREN )  COMMA  ,         
  439. C     9      9        J        T  SEMI   ;  QUOTE  '         
  440. C                      
  441.       INTEGER ALPHA(52),ALPHB(52)         
  442.       DATA ALPHA /1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,   
  443.      $    1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,1HI,1HJ,           
  444.      $    1HK,1HL,1HM,1HN,1HO,1HP,1HQ,1HR,1HS,1HT,           
  445.      $    1HU,1HV,1HW,1HX,1HY,1HZ,1H ,1H(,1H),1H;,           
  446.      $    1H:,1H+,1H-,1H*,1H/,1H\,1H=,1H.,1H,,1H',           
  447.      $    1H<,1H>/     
  448. C                      
  449. C     ALTERNATE CHARACTER SET             
  450. C                      
  451.       DATA ALPHB /1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,   
  452.      $    1Ha,1Hb,1Hc,1Hd,1He,1Hf,1Hg,1Hh,1Hi,1Hj,           
  453.      $    1Hk,1Hl,1Hm,1Hn,1Ho,1Hp,1Hq,1Hr,1Hs,1Ht,           
  454.      $    1Hu,1Hv,1Hw,1Hx,1Hy,1Hz,1H ,1H(,1H),1H;,           
  455.      $    1H|,1H+,1H-,1H*,1H/,1H$,1H=,1H.,1H,,1H",           
  456.      $    1H[,1H]/     
  457. C                      
  458.       DATA EPS/14,25,28,36/,FLOPS/15,21,24,25/               
  459.       DATA EYE/14,34,14,36/,RAND/27,10,23,13/                
  460. C                      
  461.       IF (INIT .GT. 0) GO TO 90           
  462. C                      
  463. C     RTE = UNIT NUMBER FOR TERMINAL INPUT                   
  464.       RTE = 9          
  465.       CALL FILES(RTE,BUF)                 
  466.       RIO = RTE        
  467. C                      
  468. C     WTE = UNIT NUMBER FOR TERMINAL OUTPUT                  
  469.       WTE = 9          
  470.       CALL FILES(WTE,BUF)                 
  471.       WIO = 0          
  472. C                      
  473.       IF (INIT .GE. 0) WRITE(WTE,100)     
  474.   100 FORMAT(//1X,'     < M A T L A B >'  
  475.      $  /1X,'   Version of 05/25/82')     
  476. C                      
  477. C     HIO = UNIT NUMBER FOR HELP FILE     
  478.       HIO = 11          
  479.       CALL FILES(HIO,BUF)                 
  480. C                      
  481. C     RANDOM NUMBER SEED                  
  482.       RAN(1) = 0       
  483. C                      
  484. C     INITIAL LINE LIMIT                  
  485.       LCT(2) = 25      
  486. C                      
  487.       ALFL = 52        
  488.       CASE = 0         
  489. C     CASE = 1 for file names in lower case                  
  490.       DO 20 I = 1, ALFL                   
  491.          ALFA(I) = ALPHA(I)               
  492.          ALFB(I) = ALPHB(I)               
  493.    20 CONTINUE         
  494. C                      
  495.       VSIZE = 5005     
  496.       LSIZE = 48       
  497.       PSIZE = 32       
  498.       BOT = LSIZE-3    
  499.       CALL WSET(5,0.0D0,0.0D0,STKR(VSIZE-4),STKI(VSIZE-4),1) 
  500.       CALL PUTID(IDSTK(1,LSIZE-3),EPS)    
  501.       LSTK(LSIZE-3) = VSIZE-4             
  502.       MSTK(LSIZE-3) = 1                   
  503.       NSTK(LSIZE-3) = 1                   
  504.       S = 1.0D0        
  505.    30 S = S/2.0D0      
  506.       T = 1.0D0 + S    
  507.       IF (T .GT. 1.0D0) GO TO 30          
  508.       STKR(VSIZE-4) = 2.0D0*S             
  509.       CALL PUTID(IDSTK(1,LSIZE-2),FLOPS)  
  510.       LSTK(LSIZE-2) = VSIZE-3             
  511.       MSTK(LSIZE-2) = 1                   
  512.       NSTK(LSIZE-2) = 2                   
  513.       CALL PUTID(IDSTK(1,LSIZE-1), EYE)   
  514.       LSTK(LSIZE-1) = VSIZE-1             
  515.       MSTK(LSIZE-1) = -1                  
  516.       NSTK(LSIZE-1) = -1                  
  517.       STKR(VSIZE-1) = 1.0D0               
  518.       CALL PUTID(IDSTK(1,LSIZE), RAND)    
  519.       LSTK(LSIZE) = VSIZE                 
  520.       MSTK(LSIZE) = 1                     
  521.       NSTK(LSIZE) = 1                     
  522.       FMT = 1          
  523.       FLP(1) = 0       
  524.       FLP(2) = 0       
  525.       DDT = 0          
  526.       RAN(2) = 0       
  527.       PTZ = 0          
  528.       PT = PTZ         
  529.       ERR = 0          
  530.       IF (INIT .LT. 0) RETURN             
  531. C                      
  532.    90 CALL PARSE       
  533.       IF (FUN .EQ. 1) CALL MATFN1         
  534.       IF (FUN .EQ. 2) CALL MATFN2         
  535.       IF (FUN .EQ. 3) CALL MATFN3         
  536.       IF (FUN .EQ. 4) CALL MATFN4         
  537.       IF (FUN .EQ. 5) CALL MATFN5         
  538.       IF (FUN .EQ. 6) CALL MATFN6         
  539.       IF (FUN .EQ. 21) CALL MATFN1        
  540.       IF (FUN .NE. 99) GO TO 90           
  541.       RETURN           
  542.       END
  543.                 
  544.       SUBROUTINE PARSE                    
  545.       DOUBLE PRECISION STKR(5005),STKI(5005)                 
  546.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP        
  547.       INTEGER IDS(4,32),PSTK(32),RSTK(32),PSIZE,PT,PTZ       
  548.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),HIO,RIO,WIO,RTE,WTE,FE           
  549.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)            
  550.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP          
  551.       COMMON /RECU/ IDS,PSTK,RSTK,PSIZE,PT,PTZ               
  552.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,HIO,RIO,WIO,RTE,WTE,FE                  
  553.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN  
  554.       LOGICAL EQID     
  555.       INTEGER SEMI,EQUAL,EOL,ID(4),EXCNT,LPAREN,RPAREN,COLON,PTS,ALFL           
  556.       INTEGER BLANK,COMMA,LESS,GREAT,NAME,ANS(4),ENND(4),ELSE(4),P,R            
  557.       DATA BLANK/36/,SEMI/39/,EQUAL/46/,EOL/99/,COMMA/48/,COLON/40/             
  558.       DATA LPAREN/37/,RPAREN/38/,LESS/50/,GREAT/51/,NAME/1/,ALFL/52/            
  559.       DATA ANS/10,23,28,36/,ENND/14,23,13,36/,ELSE/14,21,28,14/                 
  560. C                      
  561.    01 R = 0            
  562.       IF (ERR .GT. 0) PTZ = 0             
  563.       IF (ERR.LE.0 .AND. PT.GT.PTZ) R = RSTK(PT)             
  564.       IF (DDT .EQ. 1) WRITE(WTE,100) PT,R,PTZ,ERR            
  565.   100 FORMAT(1X,'PARSE ',4I4)             
  566.       IF (R.EQ.15) GO TO 93               
  567.       IF (R.EQ.16 .OR. R.EQ.17) GO TO 94  
  568.       SYM = EOL        
  569.       TOP = 0          
  570.       IF (RIO .NE. RTE) CALL FILES(-1*RIO,BUF)                 
  571.       RIO = RTE        
  572.       LCT(3) = 0       
  573.       LCT(4) = 2       
  574.       LPT(1) = 1       
  575.    10 IF (SYM.EQ.EOL .AND. MOD(LCT(4)/2,2).EQ.1) CALL PROMPT(LCT(4)/4)          
  576.       IF (SYM .EQ. EOL) CALL GETLIN       
  577.       ERR = 0          
  578.       PT = PTZ         
  579.    15 EXCNT = 0        
  580.       IF (DDT .EQ. 1) WRITE(WTE,115) PT,TOP                  
  581.   115 FORMAT(1X,'STATE ',2I4)             
  582.       LHS = 1          
  583.       CALL PUTID(ID,ANS)                  
  584.       CALL GETSYM      
  585.       IF (SYM.EQ.COLON .AND. CHAR.EQ.EOL) DDT = 1-DDT        
  586.       IF (SYM .EQ. COLON) CALL GETSYM     
  587.       IF (SYM.EQ.SEMI .OR. SYM.EQ.COMMA .OR. SYM.EQ.EOL) GO TO 80               
  588.       IF (SYM .EQ. NAME) GO TO 20         
  589.       IF (SYM .EQ. LESS) GO TO 40         
  590.       IF (SYM .EQ. GREAT) GO TO 45        
  591.       GO TO 50         
  592. C                      
  593. C     LHS BEGINS WITH NAME                
  594.    20 CALL COMAND(SYN)                    
  595.       IF (ERR .GT. 0) GO TO 01            
  596.       IF (FUN .EQ. 99) GO TO 95           
  597.       IF (FIN .EQ. -15) GO TO 80          
  598.       IF (FIN .LT. 0) GO TO 91            
  599.       IF (FIN .GT. 0) GO TO 70            
  600. C     IF NAME IS A FUNCTION, MUST BE RHS  
  601.       RHS = 0          
  602.       CALL FUNS(SYN)   
  603.       IF (FIN .NE. 0) GO TO 50            
  604. C     PEEK ONE CHARACTER AHEAD            
  605.       IF (CHAR.EQ.SEMI .OR. CHAR.EQ.COMMA .OR. CHAR.EQ.EOL)  
  606.      $      CALL PUTID(ID,SYN)            
  607.       IF (CHAR .EQ. EQUAL) GO TO 25       
  608.       IF (CHAR .EQ. LPAREN) GO TO 30      
  609.       GO TO 50         
  610. C                      
  611. C     LHS IS SIMPLE VARIABLE              
  612.    25 CALL PUTID(ID,SYN)                  
  613.       CALL GETSYM      
  614.       CALL GETSYM      
  615.       GO TO 50         
  616. C                      
  617. C     LHS IS NAME(...)                    
  618.    30 LPT(5) = LPT(4)                     
  619.       CALL PUTID(ID,SYN)                  
  620.       CALL GETSYM      
  621.    32 CALL GETSYM      
  622.       EXCNT = EXCNT+1                     
  623.       PT = PT+1        
  624.       CALL PUTID(IDS(1,PT), ID)           
  625.       PSTK(PT) = EXCNT                    
  626.       RSTK(PT) = 1     
  627. C     *CALL* EXPR      
  628.       GO TO 92         
  629.    35 CALL PUTID(ID,IDS(1,PT))            
  630.       EXCNT = PSTK(PT)                    
  631.       PT = PT-1        
  632.       IF (SYM .EQ. COMMA) GO TO 32        
  633.       IF (SYM .NE. RPAREN) CALL ERROR(3)  
  634.       IF (ERR .GT. 0) GO TO 01            
  635.       IF (ERR .GT. 0) RETURN              
  636.       IF (SYM .EQ. RPAREN) CALL GETSYM    
  637.       IF (SYM .EQ. EQUAL) GO TO 50        
  638. C     LHS IS REALLY RHS, FORGET SCAN JUST DONE               
  639.       TOP = TOP - EXCNT                   
  640.       LPT(4) = LPT(5)                     
  641.       CHAR = LPAREN    
  642.       SYM = NAME       
  643.       CALL PUTID(SYN,ID)                  
  644.       CALL PUTID(ID,ANS)                  
  645.       EXCNT = 0        
  646.       GO TO 50         
  647. C                      
  648. C     MULTIPLE LHS     
  649.    40 LPT(5) = LPT(4)                     
  650.       PTS = PT         
  651.       CALL GETSYM      
  652.    41 IF (SYM .NE. NAME) GO TO 43         
  653.       CALL PUTID(ID,SYN)                  
  654.       CALL GETSYM      
  655.       IF (SYM .EQ. GREAT) GO TO 42        
  656.       IF (SYM .EQ. COMMA) CALL GETSYM     
  657.       PT = PT+1        
  658.       LHS = LHS+1      
  659.       PSTK(PT) = 0     
  660.       CALL PUTID(IDS(1,PT),ID)            
  661.       GO TO 41         
  662.    42 CALL GETSYM      
  663.       IF (SYM .EQ. EQUAL) GO TO 50        
  664.    43 LPT(4) = LPT(5)                     
  665.       PT = PTS         
  666.       LHS = 1          
  667.       SYM = LESS       
  668.       CHAR = LPT(4)-1                     
  669.       CHAR = LIN(CHAR)                    
  670.       CALL PUTID(ID,ANS)                  
  671.       GO TO 50         
  672. C                      
  673. C     MACRO STRING     
  674.    45 CALL GETSYM      
  675.       IF (DDT .EQ. 1) WRITE(WTE,145) PT,TOP                  
  676.   145 FORMAT(1X,'MACRO ',2I4)             
  677.       IF (SYM.EQ.LESS .AND. CHAR.EQ.EOL) CALL ERROR(28)      
  678.       IF (ERR .GT. 0) GO TO 01            
  679.       PT = PT+1        
  680.       RSTK(PT) = 20    
  681. C     *CALL* EXPR      
  682.       GO TO 92         
  683.    46 PT = PT-1        
  684.       IF (SYM.NE.LESS .AND. SYM.NE.EOL) CALL ERROR(37)       
  685.       IF (ERR .GT. 0) GO TO 01            
  686.       IF (SYM .EQ. LESS) CALL GETSYM      
  687.       K = LPT(6)       
  688.       LIN(K+1) = LPT(1)                   
  689.       LIN(K+2) = LPT(2)                   
  690.       LIN(K+3) = LPT(6)                   
  691.       LPT(1) = K + 4   
  692. C     TRANSFER STACK TO INPUT LINE        
  693.       K = LPT(1)       
  694.       L = LSTK(TOP)    
  695.       N = MSTK(TOP)*NSTK(TOP)             
  696.       DO 48 J = 1, N   
  697.          LS = L + J-1                     
  698.          LIN(K) = IDINT(STKR(LS))         
  699.          IF (LIN(K).LT.0 .OR. LIN(K).GE.ALFL) CALL ERROR(37) 
  700.          IF (ERR .GT. 0) RETURN           
  701.          IF (K.LT.1024) K = K+1           
  702.          IF (K.EQ.1024) WRITE(WTE,47) K   
  703.    47    FORMAT(1X,'INPUT BUFFER LIMIT IS ',I4,' CHARACTERS.')                  
  704.    48 CONTINUE         
  705.       TOP = TOP-1      
  706.       LIN(K) = EOL     
  707.       LPT(6) = K       
  708.       LPT(4) = LPT(1)                     
  709.       LPT(3) = 0       
  710.       LPT(2) = 0       
  711.       LCT(1) = 0       
  712.       CHAR = BLANK     
  713.       PT = PT+1        
  714.       PSTK(PT) = LPT(1)                   
  715.       RSTK(PT) = 21    
  716. C     *CALL* PARSE     
  717.       GO TO 15         
  718.    49 PT = PT-1        
  719.       IF (DDT .EQ. 1) WRITE(WTE,149) PT,TOP                  
  720.   149 FORMAT(1X,'MACEND',2I4)             
  721.       K = LPT(1) - 4   
  722.       LPT(1) = LIN(K+1)                   
  723.       LPT(4) = LIN(K+2)                   
  724.       LPT(6) = LIN(K+3)                   
  725.       CHAR = BLANK     
  726.       CALL GETSYM      
  727.       GO TO 80         
  728. C                      
  729. C     LHS FINISHED, START RHS             
  730.    50 IF (SYM .EQ. EQUAL) CALL GETSYM     
  731.       PT = PT+1        
  732.       CALL PUTID(IDS(1,PT),ID)            
  733.       PSTK(PT) = EXCNT                    
  734.       RSTK(PT) = 2     
  735. C     *CALL* EXPR      
  736.       GO TO 92         
  737.    55 IF (SYM.EQ.SEMI .OR. SYM.EQ.COMMA .OR. SYM.EQ.EOL) GO TO 60               
  738.       IF (SYM.EQ.NAME .AND. EQID(SYN,ELSE)) GO TO 60         
  739.       IF (SYM.EQ.NAME .AND. EQID(SYN,ENND)) GO TO 60         
  740.       CALL ERROR(40)   
  741.       IF (ERR .GT. 0) GO TO 01            
  742. C                      
  743. C     STORE RESULTS    
  744.    60 RHS = PSTK(PT)   
  745.       CALL STACKP(IDS(1,PT))              
  746.       IF (ERR .GT. 0) GO TO 01            
  747.       PT = PT-1        
  748.       LHS = LHS-1      
  749.       IF (LHS .GT. 0) GO TO 60            
  750.       GO TO 70         
  751. C                      
  752. C     UPDATE AND POSSIBLY PRINT OPERATION COUNTS             
  753.    70 K = FLP(1)       
  754.       IF (K .NE. 0) STKR(VSIZE-3) = DFLOAT(K)                
  755.       STKR(VSIZE-2) = STKR(VSIZE-2) + DFLOAT(K)              
  756.       FLP(1) = 0       
  757.       IF (.NOT.(CHAR.EQ.COMMA .OR. (SYM.EQ.COMMA .AND. CHAR.EQ.EOL)))           
  758.      $       GO TO 80                     
  759.       CALL GETSYM      
  760.       I5 = 10**5       
  761.       LUNIT = WTE      
  762.    71 IF (K .EQ. 0) WRITE(LUNIT,171)      
  763.   171 FORMAT(/1X,'   no flops')           
  764.       IF (K .EQ. 1) WRITE(LUNIT,172)      
  765.   172 FORMAT(/1X,'    1 flop')            
  766.       IF (1.LT.K .AND. K.LT.100000) WRITE(LUNIT,173) K       
  767.   173 FORMAT(/1X,I5,' flops')             
  768.       IF (100000 .LE. K) WRITE(LUNIT,174) K                  
  769.   174 FORMAT(/1X,I9,' flops')             
  770.       IF (LUNIT.EQ.WIO .OR. WIO.EQ.0) GO TO 80               
  771.       LUNIT = WIO      
  772.       GO TO 71         
  773. C                      
  774. C     FINISH STATEMENT                    
  775.    80 FIN = 0          
  776.       P = 0            
  777.       R = 0            
  778.       IF (PT .GT. 0) P = PSTK(PT)         
  779.       IF (PT .GT. 0) R = RSTK(PT)         
  780.       IF (DDT .EQ. 1) WRITE(WTE,180) PT,PTZ,P,R,LPT(1)       
  781.   180 FORMAT(1X,'FINISH',5I4)             
  782.       IF (SYM.EQ.COMMA .OR. SYM.EQ.SEMI) GO TO 15            
  783.       IF (R.EQ.21 .AND. P.EQ.LPT(1)) GO TO 49                
  784.       IF (PT .GT. PTZ) GO TO 91           
  785.       GO TO 10         
  786. C                      
  787. C     SIMULATE RECURSION                  
  788.    91 CALL CLAUSE      
  789.       IF (ERR .GT. 0) GO TO 01            
  790.       IF (PT .LE. PTZ) GO TO 15           
  791.       R = RSTK(PT)     
  792.       IF (R .EQ. 21) GO TO 49             
  793.       GO TO (99,99,92,92,92,99,99,99,99,99,99,99,15,15,99,99,99,99,99),R        
  794. C                      
  795.    92 CALL EXPR        
  796.       IF (ERR .GT. 0) GO TO 01            
  797.       R = RSTK(PT)     
  798.       GO TO (35,55,91,91,91,93,93,99,99,94,94,99,99,99,99,99,99,94,94,          
  799.      $       46),R     
  800. C                      
  801.    93 CALL TERM        
  802.       IF (ERR .GT. 0) GO TO 01            
  803.       R = RSTK(PT)     
  804.       GO TO (99,99,99,99,99,92,92,94,94,99,99,99,99,99,95,99,99,99,99),R        
  805. C                      
  806.    94 CALL FACTOR      
  807.       IF (ERR .GT. 0) GO TO 01            
  808.       R = RSTK(PT)     
  809.       GO TO (99,99,99,99,99,99,99,93,93,92,92,94,99,99,99,95,95,92,92),R        
  810. C                      
  811. C     CALL MATFNS BY RETURNING TO MATLAB  
  812.    95 IF (FIN.GT.0 .AND. MSTK(TOP).LT.0) CALL ERROR(14)      
  813.       IF (ERR .GT. 0) GO TO 01            
  814.       RETURN           
  815. C                      
  816.    99 CALL ERROR(22)   
  817.       GO TO 01         
  818.       END
  819.              
  820.       SUBROUTINE PLOT(LUNIT,X,Y,N,P,K,BUF)                   
  821.       DOUBLE PRECISION X(N),Y(N),P(1)     
  822.       INTEGER BUF(79)                     
  823. C                      
  824. C     PLOT X VS. Y ON LUNIT               
  825. C     IF K IS NONZERO, THEN P(1),...,P(K) ARE EXTRA PARAMETERS                  
  826. C     BUF IS WORK SPACE                   
  827. C                      
  828.       DOUBLE PRECISION XMIN,YMIN,XMAX,YMAX,DY,DX,Y1,Y0       
  829.       INTEGER AST,BLANK,H,W               
  830.       DATA AST/1H*/,BLANK/1H /,H/20/,W/79/                   
  831. C                      
  832. C     H = HEIGHT, W = WIDTH               
  833. C                      
  834.       IF (K .GT. 0) WRITE(LUNIT,01) (P(I), I=1,K)            
  835.    01 FORMAT('Extra parameters',10f5.1)   
  836.       XMIN = X(1)      
  837.       XMAX = X(1)      
  838.       YMIN = Y(1)      
  839.       YMAX = Y(1)      
  840.       DO 10 I = 1, N   
  841.          XMIN = DMIN1(XMIN,X(I))          
  842.          XMAX = DMAX1(XMAX,X(I))          
  843.          YMIN = DMIN1(YMIN,Y(I))          
  844.          YMAX = DMAX1(YMAX,Y(I))          
  845.    10 CONTINUE         
  846.       DX = XMAX - XMIN                    
  847.       IF (DX .EQ. 0.0D0) DX = 1.0D0       
  848.       DY = YMAX - YMIN                    
  849.       WRITE(LUNIT,35)                     
  850.       DO 40 L = 1, H   
  851.          DO 20 J = 1, W                   
  852.             BUF(J) = BLANK                
  853.    20    CONTINUE      
  854.          Y1 = YMIN + (H-L+1)*DY/H         
  855.          Y0 = YMIN + (H-L)*DY/H           
  856.          JMAX = 1      
  857.          DO 30 I = 1, N                   
  858.             IF (Y(I) .GT. Y1) GO TO 30    
  859.             IF (L.NE.H .AND. Y(I).LE.Y0) GO TO 30            
  860.             J = 1 + (W-1)*(X(I) - XMIN)/DX                   
  861.             BUF(J) = AST                  
  862.             JMAX = MAX0(JMAX,J)           
  863.    30    CONTINUE      
  864.          WRITE(LUNIT,35) (BUF(J),J=1,JMAX)                   
  865.    35    FORMAT(79A1)                     
  866.    40 CONTINUE         
  867.       RETURN           
  868.       END 
  869.               
  870.       SUBROUTINE PRINT(ID,K)              
  871. C     PRIMARY OUTPUT ROUTINE              
  872.       INTEGER ID(4),K                     
  873.       DOUBLE PRECISION STKR(5005),STKI(5005)                 
  874.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP        
  875.       INTEGER ALFA(52),ALFB(52),ALFL,CASE                    
  876.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),HIO,RIO,WIO,RTE,WTE,FE           
  877.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)            
  878.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP          
  879.       COMMON /ALFS/ ALFA,ALFB,ALFL,CASE   
  880.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,HIO,RIO,WIO,RTE,WTE,FE                  
  881.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN  
  882.       DOUBLE PRECISION S,TR,TI,PR(12),PI(12),ROUND           
  883.       INTEGER FNO(11),FNL(11),SIG(12),PLUS,MINUS,BLANK,TYP,F 
  884.       DATA PLUS/41/,MINUS/42/,BLANK/36/   
  885. C     FORMAT NUMBERS AND LENGTHS          
  886.       DATA FNO /11,12,21,22,23,24,31,32,33,34,-1/            
  887.       DATA FNL /12, 6, 8, 4, 6, 3, 4, 2, 3, 1, 1/            
  888. C     FMT   1       2       3       4       5                
  889. C         SHORT   LONG   SHORT E  LONG E    Z                
  890. C     TYP   1       2       3             
  891. C         INTEGER  REAL   COMPLEX         
  892.       IF (LCT(1) .LT. 0) GO TO 99         
  893.       L = LSTK(K)      
  894.       M = MSTK(K)      
  895.       N = NSTK(K)      
  896.       MN = M*N         
  897.       TYP = 1          
  898.       S = 0.0D0        
  899.       DO 10 I = 1, MN                     
  900.         LS = L+I-1     
  901.         TR = STKR(LS)                     
  902.         TI = STKI(LS)                     
  903.         S = DMAX1(S,DABS(TR),DABS(TI))    
  904.         IF (ROUND(TR) .NE. TR) TYP = MAX0(2,TYP)             
  905.         IF (TI .NE. 0.0D0) TYP = 3        
  906.    10 CONTINUE         
  907.       IF (S .NE. 0.0D0) S = DLOG10(S)     
  908.       KS = IDINT(S)    
  909.       IF (-2 .LE. KS .AND. KS .LE. 1) KS = 0                 
  910.       IF (KS .EQ. 2 .AND. FMT .EQ. 1 .AND. TYP .EQ. 2) KS = 0                   
  911.       IF (TYP .EQ. 1 .AND. KS .LE. 2) F = 1                  
  912.       IF (TYP .EQ. 1 .AND. KS .GT. 2) F = 2                  
  913.       IF (TYP .EQ. 1 .AND. KS .GT. 9) TYP = 2                
  914.       IF (TYP .EQ. 2) F = FMT + 2         
  915.       IF (TYP .EQ. 3) F = FMT + 6         
  916.       IF (MN.EQ.1 .AND. KS.NE.0 .AND. FMT.LT.3 .AND. TYP.NE.1) F = F+2          
  917.       IF (FMT .EQ. 5) F = 11              
  918.       JINC = FNL(F)    
  919.       F = FNO(F)       
  920.       S = 1.0D0        
  921.       IF (F.EQ.21 .OR. F.EQ.22 .OR. F.EQ.31 .OR. F.EQ.32) S = 10.0D0**KS        
  922.       LS = ((N-1)/JINC+1)*M + 2           
  923.       IF (LCT(1) + LS .LE. LCT(2)) GO TO 20                  
  924.          LCT(1) = 0    
  925.          WRITE(WTE,43) LS                 
  926.          READ(RTE,44,END=19) LS           
  927. CDC..    IF (EOF(RTE).NE.0) GO TO 19      
  928.          IF (LS .EQ. ALFA(BLANK+1)) GO TO 20                 
  929.          LCT(1) = -1   
  930.          GO TO 99      
  931.    19    CALL FILES(-1*RTE,BUF)             
  932.    20 CONTINUE         
  933.       WRITE(WTE,44)    
  934.       IF (WIO .NE. 0) WRITE(WIO,44)       
  935.       CALL PRNTID(ID,-1)                  
  936.       LCT(1) = LCT(1)+2                   
  937.       LUNIT = WTE      
  938.    50 IF (S .NE. 1.0D0) WRITE(LUNIT,41) S                    
  939.       DO 80 J1 = 1, N, JINC               
  940.         J2 = MIN0(N, J1+JINC-1)           
  941.         WRITE(LUNIT,44)                   
  942.         IF (N .GT. JINC) WRITE(LUNIT,42) J1,J2               
  943.         DO 70 I = 1, M                    
  944.           JM = J2-J1+1                    
  945.           DO 60 J = 1, JM                 
  946.              LS = L+I-1+(J+J1-2)*M        
  947.              PR(J) = STKR(LS)/S           
  948.              PI(J) = DABS(STKI(LS)/S)     
  949.              SIG(J) = ALFA(PLUS+1)        
  950.              IF (STKI(LS) .LT. 0.0D0) SIG(J) = ALFA(MINUS+1) 
  951.    60     CONTINUE     
  952.           IF (F .EQ. 11) WRITE(LUNIT,11)(PR(J),J=1,JM)       
  953.           IF (F .EQ. 12) WRITE(LUNIT,12)(PR(J),J=1,JM)       
  954.           IF (F .EQ. 21) WRITE(LUNIT,21)(PR(J),J=1,JM)       
  955.           IF (F .EQ. 22) WRITE(LUNIT,22)(PR(J),J=1,JM)       
  956.           IF (F .EQ. 23) WRITE(LUNIT,23)(PR(J),J=1,JM)       
  957.           IF (F .EQ. 24) WRITE(LUNIT,24)(PR(J),J=1,JM)       
  958.           IF (F .EQ. 31) WRITE(LUNIT,31)(PR(J),SIG(J),PI(J),J=1,JM)             
  959.           IF (F .EQ. 32) WRITE(LUNIT,32)(PR(J),SIG(J),PI(J),J=1,JM)             
  960.           IF (F .EQ. 33) WRITE(LUNIT,33)(PR(J),SIG(J),PI(J),J=1,JM)             
  961.           IF (F .EQ. 34) WRITE(LUNIT,34)(PR(J),SIG(J),PI(J),J=1,JM)             
  962.           IF (F .EQ. -1) CALL FORMZ(LUNIT,STKR(LS),STKI(LS)) 
  963.           LCT(1) = LCT(1)+1               
  964.    70   CONTINUE       
  965.    80 CONTINUE         
  966.       IF (LUNIT.EQ.WIO .OR. WIO.EQ.0) GO TO 99               
  967.       LUNIT = WIO      
  968.       GO TO 50         
  969.    99 RETURN           
  970. C                      
  971.    11 FORMAT(1X,12F6.0)                   
  972.    12 FORMAT(1X,6F12.0)                   
  973.    21 FORMAT(1X,F9.4,7F10.4)              
  974.    22 FORMAT(1X,F19.15,3F20.15)           
  975.    23 FORMAT(1X,1P6D13.4)                 
  976.    24 FORMAT(1X,1P3D24.15)                
  977.    31 FORMAT(1X,4(F9.4,' ',A1,F7.4,'i'))  
  978.    32 FORMAT(1X,F19.15,A1,F18.15,'i',F20.15,A1,F18.15,'i')   
  979.    33 FORMAT(1X,3(1PD13.4,' ',A1,1PD10.4,'i'))               
  980.    34 FORMAT(1X,1PD24.15,' ',A1,1PD21.15,'i')                
  981.    41 FORMAT(/1X,' ',1PD9.1,2H *)         
  982.    42 FORMAT(1X,'    COLUMNS',I3,' THRU',I3)                 
  983.    43 FORMAT(/1X,'AT LEAST ',I5,' MORE LINES.',              
  984.      $       '  ENTER BLANK LINE TO CONTINUE OUTPUT.')       
  985.    44 FORMAT(A1)       
  986. C                      
  987.       END
  988.               
  989.       SUBROUTINE PRNTID(ID,ARGCNT)        
  990. C     PRINT VARIABLE NAMES                
  991.       INTEGER ID(4,1),ARGCNT              
  992.       INTEGER ALFA(52),ALFB(52),ALFL,CASE                    
  993.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),HIO,RIO,WIO,RTE,WTE,FE           
  994.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)            
  995.       COMMON /ALFS/ ALFA,ALFB,ALFL,CASE   
  996.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,HIO,RIO,WIO,RTE,WTE,FE                  
  997.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN  
  998.       INTEGER EQUAL    
  999.       DATA EQUAL/46/   
  1000.       J1 = 1           
  1001.    10 J2 = MIN0(J1+7,IABS(ARGCNT))        
  1002.       L = 0            
  1003.       DO 15 J = J1,J2                     
  1004.       DO 15 I = 1, 4   
  1005.       K = ID(I,J)+1    
  1006.       L = L+1          
  1007.       BUF(L) = ALFA(K)                    
  1008.    15 CONTINUE         
  1009.       IF (ARGCNT .EQ. -1) L=L+1           
  1010.       IF (ARGCNT .EQ. -1) BUF(L) = ALFA(EQUAL+1)             
  1011.       WRITE(WTE,20) (BUF(I),I=1,L)        
  1012.       IF (WIO .NE. 0) WRITE(WIO,20) (BUF(I),I=1,L)           
  1013.    20 FORMAT(1X,8(4A1,2H  ))              
  1014.       J1 = J1+8        
  1015.       IF (J1 .LE. IABS(ARGCNT)) GO TO 10  
  1016.       RETURN           
  1017.       END
  1018.              
  1019.       SUBROUTINE PROMPT(PAUSE)            
  1020.       INTEGER PAUSE    
  1021. C                      
  1022. C     ISSUE MATLAB PROMPT WITH OPTIONAL PAUSE                
  1023. C                      
  1024.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),HIO,RIO,WIO,RTE,WTE,FE           
  1025.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,HIO,RIO,WIO,RTE,WTE,FE                  
  1026.       WRITE(WTE,10)    
  1027.       IF (WIO .NE. 0) WRITE(WIO,10)       
  1028.    10 FORMAT(/1X,'<>',$)                    
  1029.       IF (PAUSE .EQ. 1) READ(RTE,20) DUMMY                   
  1030.    20 FORMAT(A1)       
  1031.       RETURN           
  1032.       END 
  1033.               
  1034.       DOUBLE PRECISION FUNCTION PYTHAG(A,B)                  
  1035.       DOUBLE PRECISION A,B                
  1036.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),HIO,RIO,WIO,RTE,WTE,FE           
  1037.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,HIO,RIO,WIO,RTE,WTE,FE                  
  1038.       DOUBLE PRECISION P,Q,R,S,T          
  1039.       P = DMAX1(DABS(A),DABS(B))          
  1040.       Q = DMIN1(DABS(A),DABS(B))          
  1041.       IF (Q .EQ. 0.0D0) GO TO 20          
  1042.       IF (DDT .EQ. 25) WRITE(WTE,1)       
  1043.       IF (DDT .EQ. 25) WRITE(WTE,2) P,Q   
  1044.     1 FORMAT(1X,'PYTHAG',1P2D23.15)       
  1045.     2 FORMAT(1X,1P2D23.15)                
  1046.    10 R = (Q/P)**2     
  1047.       T = 4.0D0 + R    
  1048.       IF (T .EQ. 4.0D0) GO TO 20          
  1049.       S = R/T          
  1050.       P = P + 2.0D0*P*S                   
  1051.       Q = Q*S          
  1052.       IF (DDT .EQ. 25) WRITE(WTE,2) P,Q   
  1053.       GO TO 10         
  1054.    20 PYTHAG = P       
  1055.       RETURN           
  1056.       END
  1057.               
  1058.       SUBROUTINE RAT(X,LEN,MAXD,A,B,D)    
  1059.       INTEGER LEN,MAXD                    
  1060.       DOUBLE PRECISION X,A,B,D(LEN)       
  1061. C                      
  1062. C     A/B = CONTINUED FRACTION APPROXIMATION TO X            
  1063. C           USING  LEN  TERMS EACH LESS THAN MAXD            
  1064. C                      
  1065.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),HIO,RIO,WIO,RTE,WTE,FE           
  1066.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,HIO,RIO,WIO,RTE,WTE,FE                  
  1067.       DOUBLE PRECISION S,T,Z,ROUND        
  1068.       Z = X            
  1069.       DO 10 I = 1, LEN                    
  1070.          K = I         
  1071.          D(K) = ROUND(Z)                  
  1072.          Z = Z - D(K)                     
  1073.          IF (DABS(Z)*DFLOAT(MAXD) .LE. 1.0D0) GO TO 20       
  1074.          Z = 1.0D0/Z   
  1075.    10 CONTINUE         
  1076.    20 T = D(K)         
  1077.       S = 1.0D0        
  1078.       IF (K .LT. 2) GO TO 40              
  1079.       DO 30 IB = 2, K                     
  1080.          I = K+1-IB    
  1081.          Z = T         
  1082.          T = D(I)*T + S                   
  1083.          S = Z         
  1084.    30 CONTINUE         
  1085.    40 IF (S .LT. 0.0D0) T = -T            
  1086.       IF (S .LT. 0.0D0) S = -S            
  1087.       IF (DDT .EQ. 27) WRITE(WTE,50) X,T,S,(D(I),I=1,K)      
  1088.    50 FORMAT(/1X,1PD23.15,0PF8.0,' /',F8.0,4X,6F5.0/(1X,45X,6F5.0))             
  1089.       A = T            
  1090.       B = S            
  1091.       RETURN           
  1092.       END              
  1093.               
  1094.       SUBROUTINE SAVLOD(LUNIT,ID,M,N,IMG,JOB,XREAL,XIMAG)    
  1095.       INTEGER LUNIT,ID(4),M,N,IMG,JOB     
  1096.       DOUBLE PRECISION XREAL(1),XIMAG(1)  
  1097. C                      
  1098. C     IMPLEMENT SAVE AND LOAD             
  1099. C     LUNIT = LOGICAL UNIT NUMBER         
  1100. C     ID = NAME, FORMAT 4A1               
  1101. C     M, N = DIMENSIONS                   
  1102. C     IMG = NONZERO IF XIMAG IS NONZERO   
  1103. C     JOB = 0     FOR SAVE                
  1104. C         = SPACE AVAILABLE FOR LOAD      
  1105. C     XREAL, XIMAG = REAL AND OPTIONAL IMAGINARY PARTS       
  1106. C                      
  1107. C     SYSTEM DEPENDENT FORMATS            
  1108.   101 FORMAT(4A1,3I4)                     
  1109.   102 FORMAT(4Z18)     
  1110. C                      
  1111.       IF (JOB .GT. 0) GO TO 20            
  1112. C                      
  1113. C     SAVE             
  1114.    10 WRITE(LUNIT,101) ID,M,N,IMG         
  1115.       DO 15 J = 1, N   
  1116.          K = (J-1)*M+1                    
  1117.          L = J*M       
  1118.          WRITE(LUNIT,102) (XREAL(I),I=K,L)                   
  1119.          IF (IMG .NE. 0) WRITE(LUNIT,102) (XIMAG(I),I=K,L)   
  1120.    15 CONTINUE         
  1121.       RETURN           
  1122. C                      
  1123. C     LOAD             
  1124.    20 READ(LUNIT,101,END=30) ID,M,N,IMG   
  1125.       IF (M*N .GT. JOB) GO TO 30          
  1126.       DO 25 J = 1, N   
  1127.          K = (J-1)*M+1                    
  1128.          L = J*M       
  1129.          READ(LUNIT,102,END=30) (XREAL(I),I=K,L)             
  1130.          IF (IMG .NE. 0) READ(LUNIT,102,END=30) (XIMAG(I),I=K,L)                
  1131.    25 CONTINUE         
  1132.       RETURN           
  1133. C                      
  1134. C     END OF FILE      
  1135.    30 M = 0            
  1136.       N = 0            
  1137.       RETURN           
  1138.       END 
  1139.               
  1140.       SUBROUTINE STACK1(OP)               
  1141.       INTEGER OP       
  1142. C                      
  1143. C     UNARY OPERATIONS                    
  1144. C                      
  1145.       DOUBLE PRECISION STKR(5005),STKI(5005)                 
  1146.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP        
  1147.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),HIO,RIO,WIO,RTE,WTE,FE           
  1148.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)            
  1149.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP          
  1150.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,HIO,RIO,WIO,RTE,WTE,FE                  
  1151.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN  
  1152.       INTEGER QUOTE    
  1153.       DATA QUOTE/49/   
  1154.       IF (DDT .EQ. 1) WRITE(WTE,100) OP   
  1155.   100 FORMAT(1X,'STACK1',I4)              
  1156.       L = LSTK(TOP)    
  1157.       M = MSTK(TOP)    
  1158.       N = NSTK(TOP)    
  1159.       MN = M*N         
  1160.       IF (MN .EQ. 0) GO TO 99             
  1161.       IF (OP .EQ. QUOTE) GO TO 30         
  1162. C                      
  1163. C     UNARY MINUS      
  1164.       CALL WRSCAL(MN,-1.0D0,STKR(L),STKI(L),1)               
  1165.       GO TO 99         
  1166. C                      
  1167. C     TRANSPOSE        
  1168.    30 LL = L + MN      
  1169.       ERR = LL+MN - LSTK(BOT)             
  1170.       IF (ERR .GT. 0) CALL ERROR(17)      
  1171.       IF (ERR .GT. 0) RETURN              
  1172.       CALL WCOPY(MN,STKR(L),STKI(L),1,STKR(LL),STKI(LL),1)   
  1173.       M = NSTK(TOP)    
  1174.       N = MSTK(TOP)    
  1175.       MSTK(TOP) = M    
  1176.       NSTK(TOP) = N    
  1177.       DO 50 I = 1, M   
  1178.       DO 50 J = 1, N   
  1179.         LS = L+MN+(J-1)+(I-1)*N           
  1180.         LL = L+(I-1)+(J-1)*M              
  1181.         STKR(LL) = STKR(LS)               
  1182.         STKI(LL) = -STKI(LS)              
  1183.    50 CONTINUE         
  1184.       GO TO 99         
  1185.    99 RETURN           
  1186.       END              
  1187.       SUBROUTINE STACK2(OP)               
  1188.       INTEGER OP       
  1189. C                      
  1190. C     BINARY AND TERNARY OPERATIONS       
  1191. C                      
  1192.       DOUBLE PRECISION STKR(5005),STKI(5005)                 
  1193.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP        
  1194.       INTEGER IDS(4,32),PSTK(32),RSTK(32),PSIZE,PT,PTZ       
  1195.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),HIO,RIO,WIO,RTE,WTE,FE           
  1196.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)            
  1197.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP          
  1198.       COMMON /RECU/ IDS,PSTK,RSTK,PSIZE,PT,PTZ               
  1199.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,HIO,RIO,WIO,RTE,WTE,FE                  
  1200.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN  
  1201.       DOUBLE PRECISION WDOTUR,WDOTUI      
  1202.       DOUBLE PRECISION SR,SI,E1,ST,E2,FLOP                   
  1203.       INTEGER PLUS,MINUS,STAR,DSTAR,SLASH,BSLASH,DOT,COLON   
  1204.       DATA PLUS/41/,MINUS/42/,STAR/43/,DSTAR/54/,SLASH/44/   
  1205.       DATA BSLASH/45/,DOT/47/,COLON/40/   
  1206. C                      
  1207.       IF (DDT .EQ. 1) WRITE(WTE,100) OP   
  1208.   100 FORMAT(1X,'STACK2',I4)              
  1209.       L2 = LSTK(TOP)   
  1210.       M2 = MSTK(TOP)   
  1211.       N2 = NSTK(TOP)   
  1212.       TOP = TOP-1      
  1213.       L = LSTK(TOP)    
  1214.       M = MSTK(TOP)    
  1215.       N = NSTK(TOP)    
  1216.       FUN = 0          
  1217.       IF (OP .EQ. PLUS) GO TO 01          
  1218.       IF (OP .EQ. MINUS) GO TO 03         
  1219.       IF (OP .EQ. STAR) GO TO 05          
  1220.       IF (OP .EQ. DSTAR) GO TO 30         
  1221.       IF (OP .EQ. SLASH) GO TO 20         
  1222.       IF (OP .EQ. BSLASH) GO TO 25        
  1223.       IF (OP .EQ. COLON) GO TO 60         
  1224.       IF (OP .GT. 2*DOT) GO TO 80         
  1225. SHAR_EOF
  1226. #    End of shell archive
  1227. exit 0
  1228. -- 
  1229. Bob Page, U of Lowell CS Dept.  page@swan.ulowell.edu  ulowell!page
  1230. Have five nice days.
  1231.